Journal of Chemical Theory and Computation

Biologically functional RNAs operate near marginal stability, and their rugged free-energy landscapes and profound structural dynamics - typically not captured by structural biology experiments - play decisive roles. Atomistic molecular dynamics (MD) simulations provide a unique means to characterize these features. However, the applicability of atomistic MD is currently limited by accessible simulation timescales and, most importantly, by force-field (FF) accuracy. Folding free energies ({Delta}G{degrees}fold) of small RNA motifs represent well-defined targets for quantitative benchmarking of RNA FFs. In practice, however, obtaining thermodynamic estimates that are sufficiently robust for direct comparison with experimental data remains highly challenging, even for small RNA systems, and many published studies rely on sampling that is not fully converged. Here, we systematically assess the performance of widely used advanced enhanced-sampling techniques using the 8-mer r(gcGAGAgc) tetraloop as a representative benchmark system. We test temperature replica exchange (T-REMD), two solute-tempering variants of replica exchange (REST2 and REHT), as well as well-tempered metadynamics and on-the-fly probability enhanced sampling combined with solute tempering (ST-MetaD and ST-OPES). Among the tested approaches, T-REMD proves to be the most robust, yielding reproducible folding equilibria and consistent estimates of {Delta}G{degrees}fold after approximately 20 s of simulation time, independent of the initial folded or unfolded conformational ensemble. Our results provide practical guidelines for selecting sampling protocols suitable for quantitative RNA benchmarks and lay the foundation for systematic validation and future refinement of RNA FFs.

The flexibility of DNA and RNA is known to play a central role in numerous biological processes, including chromatin organization and gene regulation. While a wide range of computational approaches have been developed to investigate the conformational dynamics and flexibility of proteins, analogous methods for nucleic acids remain comparatively underexplored. Elastic Network Models (ENMs) - coarse-grained mechanical representations in which macromolecules are modeled as networks of nodes connected by elastic springs - have been successfully applied to proteins, often allowing to capture experimentally observed conformational changes through a small number of harmonic normal modes. Building on a previously validated three-bead ENM for RNA, here we introduce edENM, an essential dynamics-refined ENM for DNA, RNA, and protein-nucleic acid complexes, parametrized using a diverse set of Molecular Dynamics simulations. The vibrational modes of the new edENM show good agreement with NMR data and experimental ensembles, while avoiding the unrealistic and localized deformability of previous ENM parametrizations. Additionally, we integrated this new edENM into eBDIMS, a Brownian Dynamics-based framework that enables the simulation of large-scale and anharmonic conformational transitions in protein assemblies. In this way, we are now able to explore functional motions in large protein-nucleic acid complexes such as chromatin subunits and ribosomes.

Accurate ion force fields are essential for molecular dynamics simulations of biomolecular systems, particularly in combination with modern water models such as OPC. While OPC water improves the description of bulk water and biomolecules, the transferability of existing ion force fields to this model remains an open question. Here, we systematically assess the transferability of monovalent and divalent ion force field parameters (Li+, Na+, K+, Cs+, Mg2+,Ca2+, Sr2+, Ba2+, Cl- and Br-) to OPC water by comparing single-ion and ion-pairing properties with experimental data. Our analysis reveals that no single literature parameter set provides accurate results for all ions when directly transferred to OPC water. We hence introduce the MS/G-LB(OPC) force field, which combines Mamatkulov-Schwierz-Grotz cation parameters with Loche-Bonthuis anion parameters. MS/G-LB(OPC) reproduces hydration free energies, first-shell structural properties and activity derivatives at low salt concentrations. Our results demonstrate that transferring ion parameters to OPC can lead to significant and ion-specific deviations from experimental data, making careful validation essential. At the same time, the systematic transfer and combination of ion parameters from existing force fields can provide a practical and computationally efficient alternative to full reparameterization. MS/G-LB(OPC) is available at https://git.rz.uni-augsburg.de/cbio-gitpub/opc-ion-force-fields.

Understanding RNA structural dynamics is essential for elucidating its biological functions, and molecular dynamics (MD) simulations provide an important atomistic complement to experimental approaches. However, the predictive power of MD is fundamentally limited by the accuracy of the underlying empirical Force Fields (FFs), particularly in capturing the delicate balance of non-bonded interactions. Here, we present a systematic reparameterization strategy that replaces the external gHBfix19 hydrogen-bond (H-bond) correction potential with an equivalent set of NBfix Lennard-Jones modifications within a state-of-the-art RNA FF. Using a quantitatively converged temperature replica-exchange MD ensemble of the GAGA tetraloop, we employed a reweighting-based optimization protocol to derive NBfix parameters that reproduce the thermodynamic effects of the original gHBfix19 terms. Sequential optimization of individual gHBfix19 components proved essential to ensure stable and transferable parameter refinement. The resulting fully reformulated NBfix-based variant, termed OL3CP-NBfix19, was validated on a representative set of RNA motifs, including tetranucleotides, A-form duplexes, and tetraloops. Across all tested systems, its performance is comparable to that of the reference gHBfix19 FF. By embedding the H-bond corrections directly into the standard non-bonded framework, the NBfix formulation eliminates external biasing potentials, simplifies practical deployment, and reduces computational overhead. Beyond this specific reparameterization, our results demonstrate a practical workflow for translating targeted H-bond corrections into native FF terms for efficient biomolecular simulations.

Accessing deeply buried binding sites remains a major challenge in structure-based drug discovery, where accurate description of both protein dynamics and ligand binding pathways is required. Funnel metadynamics enables simulation of complete binding processes but is computationally demanding at the all-atom resolution. By adopting the Martini 3 force field, coarse-grained funnel metadynamics (CG-FMD) substantially reduces computational requirements while retaining enhanced sampling capabilities. In this work, we assess the capability of CG-FMD to model ligand recognition at the deeply buried colchicinoids site of the tubulin {beta}-heterodimer, a multisite protein of strategic importance. We investigated the binding of colchicine, podophyllotoxin and combretastatin-A4, recovering free energy profiles with improved statistical convergence compared to AA-FMD and comparable to experimental references. In particular CG-FMD binding free energies present mean absolute errors between 3 and 10 kJ mol-1. These results propose CG-FMD as an efficient, physics-based framework for probing ligand binding to challenging sites.

The unfolding or melting temperature (TM) is a central quantity to characterize the stability of proteins and other biopolymers. The accurate prediction of protein melting temperatures by molecular mechanics force field simulations is highly desirable for many biophysical and biotechnological applications. Since the time scales for protein (un-)folding are hardly accessible in conventional MD (cMD) simulations, enhanced sampling techniques such as Temperature Replica Exchange Molecular Dynamics (TREMD) are typically employed. However, TREMD simulations are computationally very demanding especially if large temperature ranges need to be covered. Additionally, if the TM is initially unknown, setting up TREMD simulations is often challenging. To find the optimal initial conditions for such simulations, we describe their performance based on a theoretical model, which we validate on a minimalistic Markov Chain Monte Carlo (MCMC) simulation setup. In an effort to reduce the computational demand, we have investigated the possibility to use small sets of TREMD temperature ladders placed iteratively in the vicinity of a TM estimate. Different TREMD setups were extensively tested on the fast-folding protein Chignolin. We found that appropriate starting conformations lead to significantly faster convergence. Furthermore, we found that, in practice, combining multiple small temperature ladders can be advantageous in comparison to one single temperature ladder. Based on our findings, we formulate practical recommendations on how to set up TREMD for protein melting with optimal efficiency.

Molecular dynamics (MD) simulations are a powerful tool for investigating biomolecular dynamics underlying biological functions. However, the accessible spatiotemporal scales of conventional all-atom simulations remain limited by high computational costs. Coarse-graining reduces these costs by decreasing the number of interaction sites and enabling longer timesteps. In extreme cases, proteins are represented as single spherical particles; while such approximations facilitate cellular-scale simulations, they often sacrifice essential structural information, such as molecular shape and interaction anisotropy. Here, we present CGRig, a rigid-body protein model with residue-level interaction sites designed for long-time, large-scale simulations. In CGRig, each protein is treated as a single rigid-body embedding residue-level interaction sites. Its translational and rotational motions are described by the overdamped Langevin equation incorporating a shape-dependent friction matrix. Intermolecular interactions are calculated using G[o]-like native contact potentials, Debye-Huckel electrostatics, and volume exclusion. We validated that CGRig accurately reproduces the translational and rotational diffusion coefficients expected from the friction matrix for an isolated protein. For dimeric systems, the model successfully maintained native complex structures. Furthermore, two initially separated proteins converged into the correct complex with an association rate consistent with all-atom simulations. Notably, CGRig achieved a simulation performance exceeding 17 s/day for a 1,024-molecule system. These results demonstrate that CGRig provides an efficient framework for simulating protein assembly while retaining residue-level interaction specificity, making it a valuable tool for investigating large-scale biomolecular self-assembly.

Thermoresponsive phase transitions of intrinsically disordered proteins (IDPs), including both upper critical solution temperature (UCST) and lower critical solution temperature (LCST) transitions, are widely observed in natural and synthetic sequences. However, most existing coarse-grained (CG) models employ temperature-independent interactions and fail to capture solvation-driven LCST behavior. To address this, we introduce TEA (Temperature-dependent Energetics derived from hydrAtion free energies), a physics-based framework that augments CG models with temperature-dependent interactions derived from hydration thermodynamics. Leveraging extensive all-atom molecular dynamics simulations, we demonstrate that inter-residue interaction strengths quantified by excess free energies correlate linearly with hydration free energies. This correlation, combined with a validated combination rule for heterotypic interactions, allows the derivation of temperature-dependent potentials across continuous temperature space. We map these atomistic energetics to CG potentials via a single global scaling parameter, ensuring minimal overfitting and high transferability. The TEA-augmented CG models robustly distinguish UCST- and LCST-type sequences, successfully identify experimentally reported outliers, and accurately reproduce LCST-type single-chain compaction trends and phase diagrams of multiple disordered proteins. Collectively, our work provides a transferable and physically interpretable framework that bridges atomistic hydration thermodynamics and phase behavior of IDPs, enabling the simulation of thermoresponsive sequences with minimal phenomenological fitting.

Molecular conformations play a critical role in determining molecular properties such as membrane permeability, binding affinity, and ultimately therapeutic efficacy. Both experimental and computational approaches can characterize conformations within local energy minimum, and calculations of conformational free energy provide insight into why certain conformations are thermodynamically preferred over others. However, focusing solely on sampled conformations provides only a static view of the conformational landscape which may not fully illustrate why molecules result in such a conformational ensemble. Conformational transition between different conformations further explains how and why of the conformations which further inform molecule design. Here, we introduce Internal Coordinate Net (ICoN) version 1 (v1), a deep learning model trained in MD simulation data to learn the underlying physics that governed cyclic peptide conformational dynamics. ICoN-v1 enables the identification of transient conformations and all torsion rotations between local energy minima. By following the minimum-energy pathway (MEP) in the models latent space, ICoN-v1 efficiently generates fully atomistic transition pathways that capture detailed backbone and side-chain interactions governed by concerted torsional rotations. Notably, ICoN-v1 produces smooth transition pathways that are absent from the training data, demonstrating strong generalization beyond the sampled MD conformations. Analysis of the resulting concerted torsional motions and transient states highlights key residues involved at different stages of the transition, providing mechanistic insight that can inform cyclic peptide design and drug discovery.

RNA is a flexible biopolymer that adopts diverse conformations while forming structural motifs essential for its function. Classical RNA force fields often show limited transferability and inefficient sampling of transitions between stable states, particularly in moderately large RNA. To address these limitations, quantum-informed machine learning (ML) potentials have recently emerged as a promising alternative, offering improved accuracy and transferability relative to classical force fields. Here, we assess ML potentials for exploring RNA conformations using the adenine-adenine dinucleoside monophosphate (ApA) dimer, a fundamental RNA building block. We generated an extensive quantum-mechanical (QM) dataset for ApA conformations obtained from temperature replica exchange molecular dynamics (TREMD) simulations. Despite its small size, the ApA dimer exhibits six conformations in which quantum effects and solvent-mediated interactions play crucial roles. Using this dataset, we parameterized ML potentials based on the equivariant MACE architecture and informed by both ab-initio and semi-empirical data. The resulting potentials reproduce key conformational features of the ApA system, including base stacking, sugar puckering, and backbone flexibility, and provide broader coverage of structural transitions than the general-purpose SO3LR and MACE-OFF24 models. These findings highlight the importance of quantum-accurate RNA force fields towards the structural and energetic characterization of RNA complexes.

Targeted protein degradation is an emerging approach that utilizes cellular degradation pathways to inhibit a target protein. Small molecules such as molecular glues or PROTACs can be used to mediate the formation of a ternary complex with an E3 ligase and the target protein, which can dramatically enhance the degradation process. This approach is promising for cancer therapy, where degradation of oncogenic proteins can lead to cancer cell toxicity. To design new molecular glues, it is important to develop methods that predict how well a given molecule stabilizes a protein-protein interaction. However, conventional molecular dynamics simulations face challenges in capturing the long-timescale binding and unbinding events that would be used to evaluate this stabilization. In this study, we developed a strategy that allows us to evaluate the stability of protein-protein interactions in the presence of a glue molecule using weighted ensemble simulations in combination with weakened protein-protein interactions. Using this strategy, we generated unbinding trajectories of the DCAF15-RBM39 system with small molecules E7820, Indisulam, and several other Indisulam analogs. We were able to observe distinctly different behaviors between systems with different glues, which was in agreement with their reported EC50 values. We believe this approach could aid drug discovery efforts by expanding the set of druggable targets and improving the success rate of molecular glue development.

Hyperpolarization-activated cyclic nucleotide-gated (HCN) channels are are a family of voltage-gated, cyclic-nucleotide modulated Na+/K+ channels that regulate spontaneous rhythmic electrical activity in both the heart and the brain. Understanding differences in the responsiveness to cyclic adenosine monophosphate (cAMP) modulation between HCN isoforms would offer insight into the specific binding interactions that drive channel activation. Using all-atom molecular dynamics (MD) simulations and the free-energy perturbation (FEP) approach, we determined the absolute binding free energy of cAMP to the the cyclicnucleotide-binding domain (CNBD) of HCN isoforms 1-4. By studying the free-energy of ligand binding to the various isoforms of HCN, our study advances the understanding of HCN channel activation and modulation mechanisms. Overall, our work offers insight into explaining differences in channel sensitivity across the isoforms of HCN.

The catalytic mechanism of the hairpin ribozyme has remained controversial for more than two decades, with different experimental approaches often supporting distinct mechanistic interpretations. In this work, we investigate the conformational landscape of the active site along several proposed reaction pathways using all-atom molecular dynamics simulations in explicit solvent combined with enhanced sampling techniques. Specifically, we employ Hamiltonian replica exchange simulations to extensively explore active-site conformations without relying on predefined collective variables, enabling a broad characterization of the structural ensembles associated with multiple protonation states along three candidate reaction mechanisms. Our simulations suggest that a dianionic general acid/general base pathway involving direct participation of A38 and G8 is unlikely to proceed through well-defined intermediates with catalytically competent geometries. In particular, states associated with G8 deprotonation and subsequent O2 deprotonation exhibit strongly distorted active-site arrangements that appear poorly suited for reaction progression. Although highly synchronous protontransfer steps cannot be excluded, the required deprotonation of G8 remains difficult to reconcile with neutral pH conditions. In contrast, monoanionic pathways in which the non-bridging oxygens of the scissile phosphate act as transient proton relays produce intermediates that sample geometries favorable for the nucleophilic addition and leaving-group elimination steps of the reaction. These mechanisms do not require direct catalytic involvement of G8 while remaining compatible with potential acid catalysis by protonated A38+. Our results provide a unified conformational perspective on competing mechanistic scenarios. The ensembles generated here offer a foundation for future QM/MM and ML/MM calculations aimed at quantitatively resolving the free-energy landscapes governing hairpin ribozyme catalysis. Finally, the present strategy could easily be applied to other biomolecular systems with high conformational plasticity, including other ribozymes.

The MArtini Database (MAD) web server (https://mad.ens-lyon.fr/) has been updated with new tools, models, and capabilities to support a broader range of molecular systems for the Martini coarse-grained (CG) force field. The most notable addition is the MAD:Polymer Builder, enabling the automated construction of CG polymer models with varying architectures and complexities. The server now incorporates the latest developments in Martini 3, providing enhanced control within the MAD:Molecule Builder over the conversion of all-atom structures into CG representations, including the implementation of G[o]Martini 3 for protein complexes and water-protein interaction biases. The MAD:Polymer Builder and MAD:Molecule Builder are both adapted to work with intrinsically disordered proteins and domains. Substantial progress has been made in expanding the MAD:Database, making a growing library of Martini-ready compounds readily accessible across the entire MAD ecosystem. These advances position MAD as a comprehensive and evolving platform for the preparation of diverse systems in CG molecular simulations.

Artificial intelligence (AI) models have advanced rapidly, driving breakthroughs in protein structure prediction, functional annotation, and conformational exploration. Among these, molecular dynamics (MD)-inspired generative models such as AlphaFlow and BioEmu show strong potential for capturing conformational ensembles. In this study, we benchmark these models alongside physics-based MD simulations to evaluate their ability to detect cryptic pockets in proteins. Identifying such transient pockets remains a vital goal in drug discovery, as they can offer new avenues for targeting proteins traditionally challenging to modulate. We also assess two specialized residue-level predictors, PocketMiner and CryptoBank. Using the interferon inhibitory domain of Zaire Ebola VP35 (VP35), TEM-1 {beta}-lactamse with the M182T substitution (TEM {beta}-lactamase), and their mutants, we test whether each method can detect pockets and capture the effects of point mutations known to enhance or suppress pocket formation. All methods successfully identify pockets in VP35 and distinguish between opening and closing mutants. However, in TEM, where pocket opening is subtle, the methods perform inconsistently. These results highlight the promise of AI-based and simulation-based strategies in cryptic pocket discovery while pointing to the need for further improvements to achieve robust, system-independent predictions.

A variant of the U1A protein containing four substitutions to ionizable residues was generated serendipitously due to a miscommunication. Biophysical measurements show that this variant has at least twice as much helical structure as the wild-type U1A and is trimeric in solution, in contrast to the monomeric wild type. In sharp contrast, structures predicted by deep-learning AI tools (AlphaFold2 and RoseTTAFold2) and transformer-based tools (OmegaFold and ESMFold) are all highly similar to the wild-type U1A (backbone RMSD < 1 [A]). Even more surprising, two of the substituted ionizable residues are predicted to be fully buried in the non-polar core of the protein, an outcome that contradicts well-established physico-chemical principles, as ionizable residues are normally located on the protein surface. To explore this effect further, we generated sequences containing up to all twelve residues that make up the non-polar core of U1A. Across thousands of sequences, and depending on the AI model used, the majority of predicted structures contained fully buried ionizable residues while still maintaining the overall U1A fold. We then examined two additional proteins of comparable size, acylphosphatase and the de novo-designed TOP7 fold, and observed the same phenomenon: AI models frequently predicted structures with buried ionizable residues that nevertheless retained the parent fold. When these AI-predicted structures were subjected to short (50 ns) molecular dynamics simulations using physics-based force fields such as CHARMM or AMBER, the structures rapidly relaxed into ensembles that exposed ionizable residues. We conclude that while AI-based structure prediction tools perform extremely well on naturally occurring sequences, they do not reliably encode the physico-chemical principles governing the placement of ionizable residues. A straightforward remedy is to include a brief molecular dynamics simulation as a final validation step for AI-generated structures.

Coarse-grained (CG) modeling of DNA has been gathering steam in recent years due to the limitations imposed by current computer hardware, and deficiencies of atomistic force-fields. Specifically, CG simulations have emerged as a potent tool for exploring complex energy landscapes underlying biochemical processes, such as DNA transcription, as well as material design based on programmable self-assembly. In this chapter, we illustrate how the Three Interaction Site (TIS) model for DNA, a robust coarse-graining framework, can be used to study the folding landscape of DNA hairpins. We show that despite its simplicity, the TIS-DNA model quantitatively describes the hairpin folding thermodynamics and recapitulates many features of the kinetics, including the multiplicity of pathways. The free energy landscape exhibits single-funnel character with a distinct bias towards the folded state. It is likely that folding initiates through non-specific collapse of the DNA chain, involving multiple excursions on the energy landscape, until the opposing strands are approximately aligned. Subsequently, the loop region becomes more ordered, and after the first native-contact nucleates, the rest of the process becomes essentially downhill.

Molecular dynamics (MD) simulations have yielded important insights into ion conduction in potassium channels, but quantitative comparison with electrophysiological experiments remains challenging. Due to their high computational cost, MD simulations are typically performed at membrane potentials well above physiological values, and at only a limited number of voltages. Since current-voltage relationships are not necessarily linear, this limits direct comparison between simulations and experiments. Here, we introduce a method to estimate the current-voltage characteristics of ion channels from Markov state models (MSMs) constructed from MD simulations performed at only a few membrane potentials. Time-discrete MSMs of ion conduction are converted into continuous-time rate matrices, whose voltage dependence is modelled using a rate theory formulation with free energy barriers depending on membrane potential. This approach enables the prediction of channel currents over a wide voltage range without additional simulations. We validated the method using MD simulations of the potassium channels KcsA and MthK. In both cases, the currents predicted at low membrane potentials are in good agreement with those obtained directly from MD simulations, demonstrating the robustness and efficiency of the approach.

RNA function is intrinsically linked to its structural polymorphism, with molecules exploring the heterogeneous conformational ensembles resulting from complex energy landscapes. These landscapes arise from competing interactions, small energetic separations between microstates, and strong coupling to the environment, posing significant challenges for both experimental characterization and molecular simulation. In this chapter, we review current computational strategies that aim to explore RNA conformational ensembles in silico, with a specific focus on energy landscape-based approaches and atomistic simulations. We discuss key limitations related to sampling efficiency, force-field accuracy, and ensemble analysis, and illustrate their impact through case studies on a self-cleaving ribozyme and an H-type pseudoknot. Finally, we highlight emerging directions, including closer integration with experimental data and the growing role of machine learning, which will probably reinforce the predictive power of in silico RNA energy landscape exploration.

Protein kinases regulate signaling by recognizing short sequence motifs, and how these motifs bind influences both specificity and therapeutic strategies that target kinase pathways. Peptide-based inhibitors that engage substrate-recognition regions are attracting interest, but designing them requires an understanding of how a flexible peptide approaches and settles into the bound pose. Traditional studies have focused on the bound pose and affinities, whereas the steps that link the initial encounter with the bound pose have been explored less thoroughly because the relevant intermediates are too short-lived to capture experimentally and evolve on timescales that standard molecular dynamics cannot readily access. Here, we focused on Abl kinase and Abltide, the experimentally identified optimal substrate peptide for Abl kinase, and examined the sequence of events linking initial encounter to the bound pose using two-dimensional replica exchange (gREST/REUS), which selectively enhances flexibility in the peptide and its binding interface while also sampling progression along a distance coordinate. The resulting simulations yielded a detailed binding landscape, revealing five distinct encounter regions outside the substrate-binding site and six intermediate states that may connect the initial approach to the bound pose. Some encounter regions and intermediate states participate in the dominant binding pathways. During this process, EF/G/{beta}11 hydrophobic patch, together with G helix negative patch, plays a central role in guiding Abltide toward the substrate-binding site. These findings provide mechanistic insight into substrate recognition by protein kinases and offer a foundation for the rational design of peptide-based inhibitors.